Phase Coexistence in Driven One Dimensional Transport 
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We study a one-dimensional totally asymmetric exclusion process with random particle attach- 
ments and detachments in the bulk. The resulting dynamics leads to unexpected stationary regimes 
for large but finite systems. Such regimes are characterized by a phase coexistence of low and high 
density regions separated by domain walls. We use a mean-field approach to interpret the numerical 
results obtained by Monte-Carlo simulations and we predict the phase diagram of this non-conserved 
dynamics in the thermodynamic limit. 

PACS numbers: 02.50.Ey, 05.40.-a, 64.60.-i, 72.70.+m 



Even some of the simplest driven diffusive systems in 
one dimension show surprisingly rich and complex be- 
havior which is rather unexpected when looked at with 
experience gained from equilibrium phenomena yj. A 
particularly illuminating example are boundary-induced 
phase transitions in driven one-dimensional (ID) trans- 
port processes, such as the Totally Asymmetric Simple 
Exclusion Process (TASEP). The model, originally pro- 
posed in [2|, consists of particles hopping unidirection- 
ally with hard-core exclusion along a ID lattice. Due 
to conservation of the particle current in the bulk, the 
rates of incoming or outgoing particles at the boundaries 
drive the system to non-trivial stationary states Q . The 
resulting phase diagram shows continuous and discon- 
tinuous transitions of the average density of particles in 
the limit of large system sizes. These results were ob- 
tained first in mean-field theory and then extended when 
a complete analytical solution was presented solving ex- 
plicitely the recursion relations of the model or using a 
matrix product ansatz technique |4| . 

The TASEP is one out of many examples for driven 
systems with stationary non-equilibrium states, which 
cannot be described in terms of Boltzmann weights. This 
has to be contrasted with processes like the bulk adsorp- 
tion/desorption kinetics of particles on a lattice coupled 
to a reservoir ("Langmuir Kinetics" , LK), whose station- 
ary state is well described within standard concepts of 
equilibrium statistical mechanics. Here, particles adsorb 
at an empty site or desorb from an occupied one with 
fixed respective kinetic rates obeying detailed balance. 
The bulk density profile at equilibrium is described by 
a Langmuir isotherm, determined solely by the ratio of 
the two kinetic rates pf, as given by the Gibbs ensemble. 
Due to the presence of the particle reservoir there is no 
conservation of particles and no net particle current in 
the bulk. It is interesting to ask what can be expected 
in coupling two processes which have genuinely different 
dynamics and stationary states, like TASEP with open 
boundaries and LK. 

In this letter, we relax the constraint that the con- 
served dynamics in the bulk imposes to the TASEP by 



allowing particle attachment and detachment. We are 
interested in the limit where the kinetic rates are such 
that the incoming and outgoing fluxes of particles at the 
boundaries and in the bulk are comparable. This im- 
plies that a particle, injected at the boundary or attached 
somewhere in the bulk, remains long enough on the lat- 
tice to move a finite fraction of the total system size. New 
phenomena are expected in the regime of competition be- 
tween TASEP and LK for a large but finite system. Of 
course, the dynamics in an infinitely large system would 
be completely dominated by the bulk adsorption and des- 
orption rates. It turns out that the presence of the kinetic 
rates significantly change the picture of TASEP, produc- 
ing a completely reorganized phase diagram. We shall 
show by computer simulations and mean-field arguments 
that, in this non-conserved dynamics, one can have phase 
coexistence where low and high density phases are sepa- 
rated by stable discontinuities in the density profile. 

The model we discuss here is directly inspired by the 
unidirectional motion of many motor proteins along cy- 
toskeletal filaments 0. Motors advance along the fil- 
ament while attachment and detachment of motors be- 
tween the cytoplasm and the filament occur [3 ■ Recently, 
it has been shown that such dynamics can be relevant for 
modeling the filopod growth in cukaryotic cells produced 
by motor proteins interacting within actin filaments |3] • 
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FIG. 1: TASEP scheme with bulk attachment/detachment. 

We consider a ID lattice composed of sites i = l,...,N 
(Fig. ^| . The configurations are described in terms of oc- 
cupation numbers rij = 1 for a site occupied by a particle 
and Hi = for an empty site (vacancy) . The dynamics 
is determined by a master equation for the probabilities 
to find a particular configuration {rii}. We apply the 
following dynamical rules. For each time step, a site i 
is chosen at random. A particle at site i can jump to 
site i + 1 if unoccupied (we fix units of time by putting 
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this rate equal to unity). In the bulk i = 2, ..,N — 1, a 
particle can also leave the lattice with site-independent 
detachment rate wjj or fill the site (if empty) with a rate 
uja by attachment. At the boundaries, a particle can fill 
a vacancy with a rate a at site i=l, or a vacancy can be 
formed by removing a particle from the lattice with a rate 
(3 at site i — N. We refrain from giving explicitly the mas- 
ter equation for the probabilities. Correlations induced 
into the many-particle problem can be conveniently stud- 
ied within an operator representation in Fock space 
Then the equations of the bulk dynamics read: 

dn ■ 

-T£- = - ni) - Tii(l - n i+1 ) - rii) -uj D rii, 
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while at the boundaries one obtains: 
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By taking averages 0] one observes that in order to 
compute the time evolution of («i(i)) one needs the cor- 
responding averages of higher order correlations. In or- 
der to obtain an exact solution, elaborate techniques are 
necessary. We restrict the discussion to Monte-Carlo sim- 
ulations (MCS) and a mean-field approximation (MFA) 
which we shall apply below. 

The system exhibits a particle-hole symmetry in the 
following sense. A jump of a particle to the right corre- 
sponds to a vacancy move by one step to the left. Sim- 
ilarly, a particle entering the system at the left bound- 
ary can be interpreted as a vacancy leaving the lattice, 
and vice versa for the right boundary. Attachment and 
detachment of particles in the bulk is mapped to detach- 
ment and attachment of vacancies, respectively. 

We are interested in large system sizes (N ^> 1) and, 
eventually, in the "thermodynamic limit" N — > oo. In 
this case, the study of the competition between bulk and 
boundary dynamics needs that the kinetic rates decrease 
simultaneously with the system size. More precisely, we 
define the "reduced" rates ft a an d as Ha = ujaN and 
Qd=udN, keeping Ha, Qd, a, (3 fixed as N—>oo. Note 
that the binding constant K=ua/^d remains unchanged 
when passing to the thermodynamic limit. Moreover, for 
uja =ojd = 0, one arrives back at the TASEP respecting 
the same particle-hole symmetry described above. 

We have performed extensive computer simula- 
tions 17] to obtain the average density profile in the 
stationary state. We illustrate typical phenomena by fol- 
lowing a path in parameter space along curves with fixed 
a, f3, and K while increasing Qd = &a/K. Fig. |2fa) 
shows the density profile for three different values of the 
kinetic rates. At small kinetic rates, Q,a,^d *C a, j3, 
the average density (rii) in the bulk is practically con- 
stant and close to the low-density value predicted by the 
TASEP, (rii) — ol. Conversely, at high kinetic rates the 
bulk profile is structureless and essentially determined 
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FIG. 2: (a) Average density profile (rii) computed by MCS 
(continuous line) and average density profile pix) computed 
by numerical integration of MFA stationary state equations 
@ (dashed line) in the rescaled variable x = i/N for iV = 10 3 
with a = 0.2, /3 — 0.6, K — 3 and different kinetics rates Qd 
indicated in the graph, (b) MCS average density profile for dif- 
ferent system sizes, same a, f3, K as before, and Qd = 0.1. The 
width of the steep rise decreases with increasing system sizes 
N = 10 k with k = 2, 3, 4, 5 indicated in the graph, (c) MFA av- 
erage density profile for e = 10 -3 , same a,(3,K as before, and 
different kinetic rates Qd indicated in the graph. The hor- 
izontal dashed line for p(x) = 0.75 represents the Langmuir 
isotherm for K = 3. 



by the well-known ratio Kj (1 + K) of Langmuir equilib- 
rium density |l8j . A new feature appears for intermediate 
rates, fi/j = 0.1, precisely when bulk and boundary dy- 
namics compete. The density exhibits a non-monotonic 
structure in bulk, characterized by a region of low and 
high density connected by a steep rise. 

Fig. EJb) shows the density profile for different sys- 
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tern sizes. One observes a decrease of the width of the 
transition region as the number of sites is increased. The 
simulation suggests a discontinuity of the profile in terms 
of the rescaled variable x = i/N upon approaching the 
infinite system limit. A preliminary finite size scaling 
analysis is compatible with a rescaled transition width 
which scales as N~ u with v ~ 0.5. This is very differ- 
ent from the mean-field result, i^mfa = 1 .19]. Thus we 
have identified an intermediate parameter range where 
low and high density phases coexist separated by a sharp 
domain wall (DW). We also find that the discontinuity 
in the density seems to be stable or at least localized in a 
small region compared to the system size 20]. This has 
to be contrasted with the domain wall ("shock") found 
in the TASEP right at the phase boundary between the 
high and low density phase (a = (3 < 1/2) which is de- 
localized and moves as a random walker once it is far 
from the system boundaries 0. 

Phenomena like phase separation/coexistence have 
previously been observed in non-homogeneous systems 
with open boundaries like TASEP with isolated local- 
ized defects 0, 0] . The location of the domain walls 
are expected and found to be identical to the defect posi- 
tions. In contrast, the location of the DW in our homoge- 
nous model is self-tuned and determined by the values of 
the kinetic rates (see below). In systems with periodic 
boundary conditions (which are not the subject of this 
letter) phase separation has been found in TASEP with a 
blockage 0], quenched disorder 0], or in homogeneous 
systems with multi-species particle dynamics (see for 
a general criterion |16|). 

To rationalize all these findings we have developed a 
mean- field theory. Defining pi = (rii), the MFA con- 
sists of taking the average of Eqs. (|lallb|) and factorizing 
the two-site correlations, (n^n^i) = piPi + \. Then Eqs. 
(|f allb|) display the same form provided that the binary 
occupation number rii is replaced by the continuous vari- 
able pi with < pi < 1 . The equations are now interpreted 
as ordinary differential equations. 

To obtain an analytically tractable system of equations 
we have coarse-grained the discrete lattice with lattice 
constant e — L/N to a continuum. For fixed total length 
L = 1 and N — > 00, e — > one gets the nonlinear differential 
equation for the average profile in the stationary state, 

£ -d 2 xP + (2p - l)d xP + Q A (1 -p)-Q D p = 0, (2) 

where positions are measured by the rescaled variable 
x = i/N, < x < 1. Equations (|f b|) translate now to 
boundary conditions for the density field, p(0) = a and 
p(l) = l—(3. One observes that MFA respects the particle- 
hole symmetry mentioned above, provided that when 
p(x) 1 — > 1 — p(l — x) one interchanges a <-> /3, Ha <-»f2u. 
Due to this pro per ty we can restrict the discussion to the 
case ft a > 21]. The numerical mean-field solutions 
are included in Fig. E^a) for different values of flo- We 



find good agreement of MFA compared with MCS for the 
full range of kinetic rates in the limit of large N. 

In analogy with fluid dynamics, to describe these re- 
sults one considers an effective current density which for 
our problem reads j = — (e/2)d x p+p(l — p). Abbreviating 
the fluxes from and to the reservoir by T a — ^a(1 — p) and 
Td =QdP, Eq. (j2J can be read as a balance equation: 
d x j = J- a — To ■ Since there are two boundary conditions 
one has to be careful when discarding the second deriva- 
tive in for a small prefactor e. The average profile is 
then governed by similar physics as the Burgers' equa- 
tion in the inviscid limit |ll| . Generically one expects 
shocks (here, "domain walls" , DW) in the bulk and den- 
sity layers at the boundaries ("boundary layers"). Cross- 
ing a DW, the current j remains continuous in the limit 
e — > 0, while boundary layers form whenever the density 
associate to the bulk current does not fit the boundary 
condition. To better understand these features, we have 
explored the dependence of the density profile p{x) on Cljj 
for fixed a, (3 and K (see Fig. [21(c) ) . For small kinetic 
rates, f2i) = 10~ 3 , the profile is close to the one expected 
from TASEP, with a boundary layer bridging the bulk 
density up to p=l — /3 (not resolved in Fig. GJc))- In- 
creasing f^D the slope of the bulk density increases. For 
Hd > 0.05, MFA connects a region of low density (LD), 
i.e. p(x) < 1/2 to a high density region (HD), p(x) > 1/2, 
by a DW. Whereas the solution close to the left bound- 
ary is smooth, one finds a boundary layer at the right 
end bridging densities p — 1/2 down to p — I — (3. For 
larger flu the DW moves to the left, while the slope of 
the LD region increases and the HD profile flattens ap- 
proaching the Langmuir density value K/(K + 1). For 
> 1 the DW remains practically localized at the left 
boundary. Note that the DW location strongly depends 
on typical values of the bulk kinetic rates when they are 
comparable with the boundary rates a and j3. 

In the inviscid limit e — > the complete phase diagram 
can be obtained analytically within MFA, up to some 
treatment of the density discontinuities. Interestingly, 
the solution found is never given by either a constant 
low/high density profiles as in TASEP or the Lan gmu ir 
isotherm, but by a completely new set of solutions |22j . 

The mean-field analytical solution allows to draw the 
phase diagram and compare it to TASEP. Fig. El repre- 
sents a cut through the phase diagram for Qjj = 0.1 and 
K = 3 with a and (3 used as control parameters. One 
finds an extended LD-HD coexistence region separating 
a LD and a HD phase. At the boundaries of the coexis- 
tence region, the DW between the low and high density 
phases are located in the proximity of the open ends of 
the ID lattice. For small a the DW develops at the right 
end, x — 1 and moves to the left as a increases. At the 
phase boundary between the coexistence region and the 
HD phase, the DW is located at the left end of the lattice, 
x = 0. In both cases, when the DW enters and leaves the 
lattice, it matches with a boundary layer connecting the 
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FIG. 3: Phase diagrams obtained by the exact solution of 
the stationary mean-field equation @ in the inviscid limit 
for K — 3 and fiu =0.1. The inset shows the dependence of 
the DW amplitude on a for different values of j3. 

bulk density with the density given by the correspond- 
ing boundary condition. (In this case, boundary layers 
appear only for a > 1/2 and/or (3 > 1/2). MFA shows 
that the bulk density profile becomes independent of (3 
for [3 > 1/2. Hence, in this regime for a given a only 
the magnitude of the boundary layers changes, but not 
the profile of the bulk density. This explains the vertical 
phase boundaries of the coexistence region for /?>l/2. 

In addition to its location the DW is also character- 
ized by its height (see Fig. [3}. For (3 < 1/2 we find that 
the height discontinuously jumps to a finite value upon 
entering the coexistence region from the LD phase. This 
has to be contrasted with the case [3 > 1/2 where the 
phase transition from the LD to the coexistence phase is 
characterized by a continuous increase in the height of 
the DW (compare the dashed line in Fig.[3J). In MFA we 
find that both DW amplitude and position exhibit power 
law behavior with (a — ttc) 1 ^ 3 and (a — a c ) 2 / 3 , respec- 
tively. At the phase boundary to the HD phase the DW 
height always jumps to zero discontinuously. 

Working out the complete rich scenario for different £Ia 
and in the limit N 3> 1 needs a detailed analysis of 
the phase diagram ■ We just mention that the original 
maximal current phase of the TASEP appears for FIa = 
£Id only, where the Langmuir density is valued to 1/2. 
This can be proved numerically as well as analytically. 
Conversely, as soon as fi^ ^ flrj , MFA predicts that such 
a phase continuously disappears in favor of a HD phase 
if D,a>^d (respectively LD phase if £Ia<Qd)- 
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